---
title: "Stata Models in R"
author: "Emma Pendl-Robinson"
date: "10/22/2020"
output: 
  html_document:
    toc: true
    toc_float: 
      collapsed: false
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Stata code
Looking at the stata do code and executing the code

how to do mix effects model in Stata:

- offical Refeance guide: https://www.stata.com/manuals13/me.pdf

## Goal:

- replicate the stata code in R
- I did not replicate the graphs because if the tables are the same then the graphs should be the same


```{r, message=FALSE, warning=FALSE}
library(haven)
library(tidyverse)
library(lme4)
library(magrittr)
```


Stata:
```
use "imm.bjpols.dta"

```

R:
```{r}
df <- read_dta ("../../data/raw/imm.bjpols.dta")

# what varriables are categorical 
#fact_vars <- c('uid', 'country','gender','education2','occupation2','treat_hstat', 'treat_drk', 'treat_nat_me', 'treat_kids', 'cand', 'treat_gender', 'insample','inc3')

fact_vars <- c('cand')

df %<>% mutate_at(fact_vars, as.factor)
```

# FIRST MODEL: BASIC RESULTS: ALL TREATMENTS, NO INTERACTIONS

## ALL COUNTRIES COMBINED
Page 1 of PDF stata log
Stata:
```
set more off <- this suppresses some of the warning messages
mixed support i.treat_hstat i.treat_drk i.treat_nat_me i.treat_kids i.cand if country~=12 & country~=6 & treat_gender==0 || id: [weight=wt_country] 
```
- `set more off` - suppresses some of the warning messages  
- `mixed` - Multilevel mixed-effects linear regression

R:
- Note: my output coeffients are very close. 

Resources:
- https://stats.stackexchange.com/questions/12709/multilevel-regression-using-lmer-function-in-r-and-stata
- http://www-personal.umich.edu/~bwest/mccoach_etal_2018_mlmcompare.pdf
```{r}
#filter the data to where country does not equal 12 or 6 and treat_gender is 0
df_filtered <- df %>% filter(country !=  12 & country != 6 & treat_gender ==0)

first_all_countries <- lmer(support ~ treat_hstat + treat_drk + treat_nat_me + treat_kids + cand + (1|id),
                            data = df_filtered,
                            REML = FALSE,
                            weights = wt_country
                            )

summary(first_all_countries)
```

Stata:
```
estimates store Aa
```
- I think `estimates` is just storing the estimated values https://www.stata.com/manuals13/restimatesstore.pdf

```
margins i.treat_drk 
marginsplot, title("All Countries", size(large) position(11) ring(0) margin(2 2 2 2)) scheme(s1mono) ytitle("Level of Support") xtitle("") ///
	xscale(range(-.3 1.3)) xlab(0 "Light" 1 "Dark") ///
	yscale(range(.3 .7)) ylab(.3 .4 .5 .6 .7, angle(0)) yline(0, lstyle(dot)) ///
	plotopts( msymbol(S) mcolor(black)) ciopts(lpattern(dash)) plotregion(lstyle(none)) nodraw
graph save Aa2, replace
margins i.treat_nat_me 
marginsplot, title("All Countries", size(large) position(11) ring(0) margin(2 2 2 2)) scheme(s1mono) ytitle("Level of Support") xtitle("") ///
	xscale(range(-.3 1.3)) xlab(0 "Not ME" 1 "ME") ///
	yscale(range(.3 .7)) ylab(.3 .4 .5 .6 .7, angle(0)) yline(0, lstyle(dot)) ///
	plotopts( msymbol(S) mcolor(black)) ciopts(lpattern(dash)) plotregion(lstyle(none)) nodraw
graph save Aa3, replace

```
- The margins command estimates margins of responses for specified values of covariates and presents the results as a table https://www.stata.com/manuals13/rmargins.pdf


Notes:
- I can not replicate the margins in R. I think this is because the indepent varriables are continious or maybe because things are supposed to be on the log scale? 
```{r}
library(margins)

summary(margins(first_all_countries))

```

- https://cran.r-project.org/web/packages/margins/vignettes/Stata.html
-https://cran.r-project.org/web/packages/margins/vignettes/Introduction.html
 
- http://sia.webpopix.org/EMlme.html



## BY COUNTRY

Stata: 
```
*BY COUNTRY
xtset uid
set more off

foreach C in "AU" "CA" "CH" "DK" "ES" "FR" "JP" "KR" "NO" "UK" "US" {
	xtreg support i.treat_hstat i.treat_drk i.treat_nat_me i.treat_kids i.cand  if countryt=="`C'" & treat_gender==0, mle
	estimates store A_`C'
		margins i.treat_hstat 
	marginsplot, title("`C'", size(large) position(11) ring(0) margin(2 2 2 2)) scheme(s1mono) ytitle("Level of Support") xtitle("") ///
		xscale(range(-.3 1.3)) xlab(0 "Low Status" 1 "High Status") ///
		yscale(range(.3 .7)) ylab(.3 .4 .5 .6 .7, angle(0)) yline(0, lstyle(dot)) ///
		plotopts( msymbol(S) mcolor(black)) ciopts(lpattern(dash)) plotregion(lstyle(none)) nodraw 
	graph save Aa1_`C', replace
	margins i.treat_drk 
	marginsplot, title("`C'", size(large) position(11) ring(0) margin(2 2 2 2)) scheme(s1mono) ytitle("Level of Support") xtitle("") ///
		xscale(range(-.3 1.3)) xlab(0 "Light" 1 "Dark") ///
		yscale(range(.3 .7)) ylab(.3 .4 .5 .6 .7, angle(0)) yline(0, lstyle(dot)) ///
		plotopts( msymbol(S) mcolor(black)) ciopts(lpattern(dash)) plotregion(lstyle(none)) nodraw
	graph save Aa2_`C', replace
	margins i.treat_nat_me 
	marginsplot, title("`C'", size(large) position(11) ring(0) margin(2 2 2 2)) scheme(s1mono) ytitle("Level of Support") xtitle("") ///
		xscale(range(-.3 1.3)) xlab(0 "Not ME" 1 "ME") ///
		yscale(range(.3 .7)) ylab(.3 .4 .5 .6 .7, angle(0)) yline(0, lstyle(dot)) ///
		plotopts( msymbol(S) mcolor(black)) ciopts(lpattern(dash)) plotregion(lstyle(none)) nodraw
	graph save Aa3_`C', replace
	}
```
- `xtset` is some king of memory thing
- runing for loop to loop through all the countries running `xtreg` with `mle`, finding the `margins` and making a graph
- `xtreg` fits regression models to panel data. In particular, xtreg with the be option fits random effects models by using the between regression estimator https://www.stata.com/manuals13/xtxtreg.pdf
- we want to speficy that we want to run a full maximum likelihood estimate vs a gmm https://www.stata.com/support/faqs/statistics/xtreg-mle-versus-gmm/



R:
Notes
- I think the `plm` outputs are almost identical to the `xtreg`  
    - `_cons` in stata is `(Intercept)` in r
    - I think `cand` is the order of the vignettes
    - TODO: find out what this error message means `indexes of pseries have same length but not same content: result was assigned first operand's index`
        - This is typically shown when the index's are not unique pairs but I know that they are unique pairs https://stackoverflow.com/questions/43663594/error-in-plm-regression/43664969
- I am not sure how Valentino ran `margin` to get marginal effects for pannel data. I don't think this is an option in R and I am not sure if it is mathematically possible. I think Valentino might have predicted marginal effects for  `cand` {1, 2} separetly but I am not sure
    - TODO: read more about marginal effects and how to use them
- I am going to leave out the plots
```{r}
library("plm")

#countryt_ls <- c('AU', 'CA', 'CH','DK', 'ES', 'FR', 'JP', 'KR', 'NO', 'UK', 'US')
countryt_ls <- c('AU', 'CA', 'DK', 'FR', 'JP', 'KR', 'NO', 'UK', 'US')

for (country_name in countryt_ls) {
  print(country_name)
  df_filtered <- df %>% filter(countryt == country_name & treat_gender ==0)
  
  first_by_country <- plm(support ~ treat_hstat + treat_drk + treat_nat_me + treat_kids + cand, 
                        data = df_filtered,
                        model="random", 
                        index=c("uid","cand"))
  
  print(summary(first_by_country))
  #summary(margins(first_by_country))
}

```
- https://stackoverflow.com/questions/59923956/trying-to-reproduce-xtreg-in-stata-with-plm-in-r
- https://www.princeton.edu/~otorres/Panel101R.pdf


# SECOND MODEL: ADDING INTERACTIONS

## ALL COUNTRIES COMBINED
- page 27 in pdf version of the log

Stata:
- 

```
mixed support i.treat_hstat##i.treat_drk i.treat_nat_me i.treat_hstat##i.treat_kids i.cand if country~=12 & country~=6 & treat_gender==0 || id: [weight=wt_country] 

```
R:
Notes
- most of the coefficents look ok but the interaction effects are off
    - I think this is an issue with me misspecifying `lmer` because `plm` is working fine
```{r}
#filter the data to where country does not equal 12 or 6 and treat_gender is 0
df_filtered <- df %>% filter(country !=  12 & country != 6 & treat_gender ==0)

sec_all_countries <- lmer(support ~ treat_hstat*treat_drk + treat_nat_me + treat_hstat*treat_kids + cand + (1|id),
                            data = df_filtered,
                            REML = FALSE,
                            weights = wt_country
                            )

summary(sec_all_countries)
```

- https://stats.idre.ucla.edu/r/seminars/interactions-r/


Stata:
```
margins i.treat_hstat#treat_kids

```
R:
- I don't think I am using margins correctly
```{r}
summary(margins(sec_all_countries))
```

## BY COUNTRY

Stata:
```
foreach C in "AU" "CA" "CH" "DK" "ES" "FR" "JP" "KR" "NO" "UK" "US" {
	xtreg support i.treat_hstat##i.treat_drk i.treat_nat_me i.treat_hstat##i.treat_kids i.cand  if countryt=="`C'" & treat_gender==0, mle
	estimates store B_`C'
	margins i.treat_hstat#treat_kids 
	marginsplot, title("`C'", size(large) position(11) ring(0) margin(2 2 2 2)) scheme(s1mono) ytitle("Level of Support") xtitle("") ///
		xscale(range(-.3 1.3)) xlab(0 "Low Status" 1 "High Status") ///
		yscale(range(.3 .7)) ylab(.3 .4 .5 .6 .7, angle(0)) yline(0, lstyle(dot)) ///
		plotopts( msymbol(S) mcolor(black)) ciopts(lpattern(dash)) ///
		legend(off) ysize(3.66) xsize(3.7) plotregion(lstyle(none)) nodraw
	graph save Bb1_`C', replace
	}
```
R:
- model coeffients are the same
- TODO: write code for margines
- TODO fix issue with matrix singularity 
```{r}
#countryt_ls <- c('AU', 'CA', 'CH','DK', 'ES', 'FR', 'JP', 'KR', 'NO', 'UK', 'US')
countryt_ls <- c('AU', 'CA','DK', 'FR', 'JP', 'KR', 'NO', 'UK', 'US')

for (country_name in countryt_ls) {
  print(country_name)
  df_filtered <- df %>% filter(countryt == country_name & treat_gender ==0)
  
  sec_by_country <- plm(support ~ treat_hstat*treat_drk + treat_nat_me + treat_hstat*treat_kids + cand, 
                        data = df_filtered,
                        model="random", 
                        index=c("uid","cand"))
  
  print(summary(sec_by_country))
  #summary(margins(sec_by_country))
}
```

# THIRD MODEL: ADDING SES

## ALL COUNTRIES COMBINED
page 49 of pdf copy of Stata log

### JOB STATUS BY EDUCATION

Stata:
```
mixed support i.treat_hstat##i.treat_drk i.treat_hstat##i.treat_kids##i.education2 i.treat_nat_me i.cand if country~=12 & country~=6 & treat_gender==0 || id: [weight=wt_country] 
```
R:
- my interaction coeffiecents are close but not the same
```{r}
#filter the data to where country does not equal 12 or 6 and treat_gender is 0
df_filtered <- df %>% filter(country !=  12 & country != 6 & treat_gender ==0)

third_all_job_education <- lmer(support ~ treat_hstat*treat_drk + treat_hstat*treat_kids*education2 + treat_nat_me + cand + (1|id),
                            data = df_filtered,
                            REML = FALSE,
                            weights = wt_country
                            )

#summary(third_all_job_education)
```



Stata:
```
margins treat_hstat#education2
```
R:
- TODO: find out how to do margins
```{r}
#summary(margins(third_all_job_education))
```


### JOB STATUS BY OCCUPATION
page 51

Stata:
```
mixed support i.treat_hstat##i.treat_drk i.treat_hstat##i.treat_kids##i.occupation2 i.treat_nat_me i.cand if country~=12 & country~=6 & treat_gender==0 || id: [weight=wt_country]
```

R:
- correlation coeffecnts are a little off
```{r}
#filter the data to where country does not equal 12 or 6 and treat_gender is 0
df_filtered <- df %>% filter(country !=  12 & country != 6 & treat_gender ==0)

third_all_job_occupation <- lmer(support ~ treat_hstat*treat_drk + treat_hstat*treat_kids*occupation2 + treat_nat_me + cand + (1|id),
                            data = df_filtered,
                            REML = FALSE,
                            weights = wt_country
                            )

#summary(third_all_job_occupation)
```

### JOB STATUS BY INCOME 
page 53

Stata:
```
mixed support i.treat_hstat##i.treat_drk i.treat_hstat##i.treat_kids##i.inc2 i.treat_nat_me i.cand if country~=12 & country~=6 & treat_gender==0 || id: [weight=wt_country] 

```

R:
- correlation coeffecnts are a little off
```{r}
#filter the data to where country does not equal 12 or 6 and treat_gender is 0
df_filtered <- df %>% filter(country !=  12 & country != 6 & treat_gender ==0)

third_all_job_income <- lmer(support ~ treat_hstat*treat_drk + treat_hstat*treat_kids*inc2 + treat_nat_me + cand + (1|id),
                            data = df_filtered,
                            REML = FALSE,
                            weights = wt_country
                            )

#summary(third_all_job_income)
```

## BY COUNTRY
Page 55

### JOB STATUS BY EDUCATION

Stata:
```
foreach C in "AU" "CA" "CH" "DK" "ES" "FR" "JP" "KR" "NO" "UK" "US" {
	xtreg support i.treat_hstat##i.treat_drk i.treat_hstat##i.treat_kids##i.education2 i.treat_nat_me i.cand if countryt=="`C'" & treat_gender==0, mle
	estimates store C_`C'
	margins treat_hstat#education2 
	marginsplot, title("`C'", size(large) position(11) ring(0) margin(2 2 2 2)) scheme(s1mono) ytitle("Level of Support") xtitle("") ///
		xscale(range(-.3 1.3)) xlab(0 "Low Status" 1 "High Status") ///
		yscale(range(.3 .7)) ylab(.3 .4 .5 .6 .7, angle(0)) yline(0, lstyle(dot)) ///
		plotopts( msymbol(S) mcolor(black)) ciopts(lpattern(dash))  ///
		legend(off) plotregion(lstyle(none)) nodraw
	graph save Cc1_`C', replace
	}

```

R:
Notes:
- I get the error message `Error in solve.default(crossprod(ZBeta)) : Lapack routine dgesv: system is exactly singular: U[12,12] = 0` I think this means that there is a singularity issue (i.e. correlation is too strong or not enought varraince).
- TODO: fix the singularity issue for ES
```{r}
#countryt_ls <- c('AU', 'CA', 'CH','DK', 'ES', 'FR', 'JP', 'KR', 'NO', 'UK', 'US')
countryt_ls <- c('AU', 'CA', 'CH','DK', 'FR', 'JP', 'KR', 'NO', 'UK', 'US')

# removing rows where two vingets were not completed 
df_by_country <- df %>%
  filter(treat_gender ==0) %>%
  select(countryt, treat_gender, support, treat_hstat, treat_drk, treat_kids, education2, treat_nat_me, cand, uid) %>%
  na.omit() %>%
  group_by(uid) %>% 
  mutate(n=n()) %>% 
  filter(n == 2) %>% 
  ungroup()

third_by_country_job_education <- list() # create an empty object to store results

for (country_name in countryt_ls) {
  df_filtered <- df_by_country %>% 
    filter(countryt == country_name & treat_gender == 0) 
  
  country_plm <- plm(support ~ treat_hstat*treat_drk  + treat_hstat*treat_kids*education2 + treat_nat_me + cand, 
                        data = df_filtered,
                        model="random", 
                        index=c("uid","cand"))
  
  third_by_country_job_education[[country_name]] <- country_plm
}
```

```{r eval=FALSE}
#ES is sigularity issue Page 63
# the data frame is the correct size
df_filtered <- df %>% filter(countryt == "ES" & treat_gender ==0) %>%
  select(support, treat_hstat, treat_drk, treat_kids, education2, treat_nat_me, cand, uid) %>%
  na.omit() %>%
  group_by(uid) %>% 
  mutate(n=n()) %>% 
  filter(n == 2) %>% 
  ungroup()
  

  test <- plm(support ~ treat_hstat*treat_drk  + treat_hstat*treat_kids*education2 + treat_nat_me + cand,
                        data = df_filtered,
                        model="random", 
                        index=c("uid","cand"))
  
  summary(test)
```

```{r, results='asis'}
stargazer(third_all_job_education, third_by_country_job_education[c('AU', 'CA', 'DK', 'FR', 'JP')],
          title = "Appendix Table 3. Adding SES interactions (Education)", 
          column.labels = c("All", 'AU', 'CA', 'DK', 'FR', "JP"),
          df = FALSE,
          model.names	=TRUE, model.numbers = FALSE,
          font.size = "footnotesize", 
          column.sep.width = "1pt", 
          single.row = FALSE,
          digits = 2,
          no.space = TRUE)

stargazer(third_by_country_job_education[c('KR', 'NO', 'CH', 'UK', 'US')],
          title = "Appendix Table 3 (Continued). Adding SES interactions (Education)", 
          column.labels = c('KR', 'NO', 'CH', 'UK', 'US'),
          df = FALSE,
          model.names	=TRUE, model.numbers = FALSE,
          font.size = "footnotesize", 
          column.sep.width = "1pt", 
          single.row = FALSE,
          digits = 2,
          no.space = TRUE)
```

### JOB STATUS BY OCCUPATION
Page 78

Stata:
```
foreach C in "AU" "CA" "DK" "FR" "JP" "KR" "UK" "US" {
	xtreg support i.treat_hstat##i.treat_drk i.treat_hstat##i.treat_kids##i.occupation2 i.treat_nat_me i.cand  if countryt=="`C'" & treat_gender==0, mle
	estimates store D_`C'
	}	
```


R:
Notes:
- I get the error message `Error in solve.default(crossprod(ZBeta)) : system is computationally singular: reciprocal condition number = 1.19889e-18` I think this means that I can not invert my matrix.
- TODO: fix the singularity issue for CH, ES, and NO - I don't think I need to fix this because these countries are not in the table 4 graphs
```{r}
#countryt_ls <- c('AU', 'CA', 'CH','DK', 'ES', 'FR', 'JP', 'KR', 'NO', 'UK', 'US')
countryt_ls <- c('AU', 'CA','DK', 'FR', 'JP', 'KR', 'UK', 'US')

# removing rows where two vingets were not completed 
df_by_country <- df %>%
  filter(treat_gender ==0) %>%
  select(countryt, treat_gender, support, treat_hstat, treat_drk, treat_kids, occupation2, treat_nat_me, cand, uid) %>%
  na.omit() %>%
  group_by(uid) %>% 
  mutate(n=n()) %>% 
  filter(n == 2) %>% 
  ungroup()

third_by_country_job_occupation <- list() # create an empty object to store results

for (country_name in countryt_ls) {
  df_filtered <- df_by_country %>% filter(countryt == country_name & treat_gender ==0)
  country_plm <- plm(support ~ treat_hstat*treat_drk  + treat_hstat*treat_kids*occupation2 + treat_nat_me + cand, 
                        data = df_filtered,
                        model="random", 
                        index=c("uid","cand"))
  
  third_by_country_job_occupation[[country_name]] <- country_plm
}
```


```{r, results='asis'}
stargazer(third_all_job_occupation, third_by_country_job_occupation[c('AU', 'CA', 'DK', 'FR', 'JP')],
          title = "Appendix Table 4. Adding SES interactions (Occupation)", 
          column.labels = c("All", 'AU', 'CA', 'DK', 'FR', "JP"),
          df = FALSE,
          model.names	=TRUE, model.numbers = FALSE,
          font.size = "footnotesize", 
          column.sep.width = "1pt", 
          single.row = FALSE,
          digits = 2,
          no.space = TRUE)

stargazer(third_by_country_job_occupation[c('KR', 'UK', 'US')],
          title = "Appendix Table 4 (Continued). Adding SES interactions (Occupation)", 
          column.labels = c('KR','UK', 'US'),
          df = FALSE,
          model.names	=TRUE, model.numbers = FALSE,
          font.size = "footnotesize", 
          column.sep.width = "1pt", 
          single.row = FALSE,
          digits = 2,
          no.space = TRUE)
```

### JOB STATUS BY INCOME
page 90
Stata:
```
foreach C in "AU" "CA" "DK" "FR" "JP" "KR" "ES" "CH" "UK" "US" {
	xtreg support i.treat_hstat##i.treat_drk i.treat_hstat##i.treat_kids##i.inc2 i.treat_nat_me i.cand  if countryt=="`C'" & treat_gender==0, mle
	estimates store E_`C'
	}
```
R:
Notes:
- I get the error message `Error in solve.default(crossprod(ZBeta)) : Lapack routine dgesv: system is exactly singular: U[12,12] = 0` I think this means that there is a singularity issue (i.e. correlation is too strong or not enought varraince).
- TODO: fix the singularity issue for ES and NO - do not have to fix NO
```{r}
#countryt_ls <- c('AU', 'CA', 'CH','DK', 'ES', 'FR', 'JP', 'KR', 'NO', 'UK', 'US')
countryt_ls <- c('AU', 'CA', 'CH','DK', 'FR', 'JP', 'KR', 'UK', 'US')

# removing rows where two vingets were not completed 
df_by_country <- df %>%
  filter(treat_gender ==0) %>%
  select(countryt, treat_gender, support, treat_hstat, treat_drk, treat_kids, inc2, treat_nat_me, cand, uid) %>%
  na.omit() %>%
  group_by(uid) %>% 
  mutate(n=n()) %>% 
  filter(n == 2) %>% 
  ungroup()

third_by_country_job_income <- list() # create an empty object to store results

for (country_name in countryt_ls) {
  df_filtered <- df_by_country %>% filter(countryt == country_name & treat_gender ==0)
  country_plm <- plm(support ~ treat_hstat*treat_drk  + treat_hstat*treat_kids*inc2 + treat_nat_me + cand, 
                        data = df_filtered,
                        model="random", 
                        index=c("uid","cand"))
  
  third_by_country_job_income[[country_name]] <- country_plm
}
```

```{r, results='asis'}
stargazer(third_all_job_income, third_by_country_job_income[c('AU', 'CA', 'DK', 'FR', 'JP')],
          title = "Appendix Table 5. Adding SES interactions (Income)", 
          column.labels = c("All", 'AU', 'CA', 'DK', 'FR', "JP"),
          df = FALSE,
          model.names	=TRUE, model.numbers = FALSE,
          font.size = "footnotesize", 
          column.sep.width = "1pt", 
          single.row = FALSE,
          digits = 2,
          no.space = TRUE)

stargazer(third_by_country_job_income[c('KR', 'CH', 'UK', 'US')],
          title = "Appendix Table 5 (Continued). Adding SES interactions (Income)", 
          column.labels = c('KR', 'CH', 'UK', 'US'),
          df = FALSE,
          model.names	=TRUE, model.numbers = FALSE,
          font.size = "footnotesize", 
          column.sep.width = "1pt", 
          single.row = FALSE,
          digits = 2,
          no.space = TRUE)
```

# The Impact of Economic Concerns on Openness to Immigration
page 106

Stata:
- I think this is a straight up regression https://www.stata.com/manuals13/rregress.pdf
```
reg openness c.inctaxes##i.occupation2 c.takejobs##i.occupation2 i.country if country~=12 & country~=6 & cand==1 [weight=wt_country]
```

R:
- coefficents look the same

```{r}
#filter the data to where country does not equal 12 or 6 and vignette number (cand) ==1
df_filtered <- df %>% filter(country !=  12 & country != 6 & cand ==1)

opennes_occupation <- lm(openness ~ inctaxes*occupation2 + takejobs*occupation2 + countryt,
                            data = df_filtered,
                            weights = wt_country)
summary(opennes_occupation)
```

Stata:
```
reg openness c.inctaxes##i.inc2 c.takejobs##i.inc2 i.country if country~=12 & country~=6 & cand==1 [weight=wt_country]
```

R:
- coefficents look the same

```{r}
#filter the data to where country does not equal 12 or 6 and vignette number (cand) ==1
df_filtered <- df %>% filter(country !=  12 & country != 6 & cand ==1)

opennes_income <- lm(openness ~ inctaxes*inc2 + takejobs*inc2 + countryt,
                            data = df_filtered,
                            weights = wt_country)
summary(opennes_income)
```

	
# FOURTH MODEL: ADDING RACIAL ANIMUS
## BY COUNTRY
Page 108
Stata:
```
foreach C in "CA" "FR" "ES" "UK" "US" {
	xtreg support i.treat_hstat##i.treat_drk i.treat_nat_me i.treat_hstat##i.treat_kids i.cand animus gender age i.education2 i.inc3 if countryt=="`C'" & treat_gender==0, mle
	estimates store G_`C'	
	xtreg support i.treat_hstat##i.treat_drk i.treat_nat_me i.treat_hstat##i.treat_kids i.cand i.treat_hstat##c.animus gender age i.education2 i.inc3 if countryt=="`C'" & treat_gender==0, mle
	estimates store H_`C'
	}	
```

R:
Notes:
- TODO add estimates
- TODO: address error message: `Error in solve.default(vcov(x)[names(coefs_wo_int), names(coefs_wo_int)], : Lapack routine dgesv: system is exactly singular: U[1,1] = 0` for `ES`
```{r}
#countryt_ls <- c('CA', 'FR', 'ES', 'UK', 'US')
countryt_ls <- c('CA', 'FR', 'UK', 'US')

for (country_name in countryt_ls) {
  print(country_name)
  df_filtered <- df %>% filter(countryt == country_name & treat_gender ==0) %>% mutate(inc3 = as.factor(inc3))
  
  #anumius no interaction
  fourth_by_country <- plm(support ~ treat_hstat*treat_drk  + treat_nat_me + treat_hstat*treat_kids + cand + animus + gender + age + education2 + inc3,
                        data = df_filtered,
                        model="random", 
                        index=c("uid","cand"))
  
  print(summary(fourth_by_country))
  #summary(margins(fourth_by_country))
  
  #anumius with interaction
  fourth_by_country_w_interaction <- plm(support ~ treat_hstat*treat_drk  + treat_nat_me + treat_hstat*treat_kids + cand + treat_hstat*animus + gender + age + education2 + inc3,
                        data = df_filtered,
                        model="random", 
                        index=c("uid","cand"))
  
  print(summary(fourth_by_country))
  #summary(margins(fourth_by_country))

}
```
